!------------------------------------------------------------------------------
! TJU/Department of Mechanics, Fluid Mechanics, Code START
!------------------------------------------------------------------------------
!> @Author: Jianxin Liu shookware@tju.edu.cn
!> @Date: 2024-08-05 10:35:15
!> @LastEditTime: 2024-08-05 10:35:15
!> @LastEditors: Jianxin Liu shookware@tju.edu.cn
!> @Description:
!> @FilePath: /START_NS_3D/test/field/ex2.f90
!> @Copyright (c) 2024 by Jianxin Liu email: shookware@tju.edu.cn, All Rights Reserved.
!------------------------------------------------------------------------------

program ex_field_conservative
   use petsc_wrapper_system
   use petsc_wrapper_vector
   use petsc_wrapper_domain
   use mod_field
   use mod_eos_base
   use petsc
   implicit none

   type(type_domain) :: domain, domain1
   type(type_domain_vector) :: Q
   type(type_conservative) :: conservative
   type(type_primitive) :: primitive
   type(type_domain_vector) :: u_vec, rho_vec, v_vec, w_vec, T_vec
   integer :: i, ierr

   type(type_eos) :: eos
   real(8) :: T = 300.0d0, rho = 1.293d0, u = 200.0d0, v = 0.0d0, w = 0.0d0, R = 287.15
   real(8) :: T_, rho_, u_, v_, w_, rho_e_

   real(8) :: Ma = 6.0d0

   call petsc_initialize()
   call eos%create("ideal_gas")
   call eos%ref%set("r", R)
   call eos%ref%set("gamma", eos%gas%get_cp()/eos%gas%get_cv())
   call eos%ref%set("velocity", u)
   call eos%ref%set("temperature", T)
   call eos%ref%set("density", rho)
   call eos%ref%set("flg_uref", 1.0d0/Ma)
   T_ = T/eos%ref%T
   rho_ = rho/eos%ref%rho
   u_ = u/eos%ref%velocity
   v_ = v/eos%ref%velocity
   w_ = w/eos%ref%velocity

   print *, rho_, u_, v_, w_, T_

   print *, eos%ref
   rho_e_ = 1.0d0/eos%ref%gamma*1.0d0/Ma**2*eos%gas%get_cv()/eos%ref%R*rho_*T_ + 0.5d0*rho_*(u_**2 + v_**2 + w_**2)
   call domain%create([5], dof=5, stencil_width=1)
   call Q%create(domain)
   call Q%set_value_blocked([(rho_, rho_*u_, rho_*v_, rho_*w_, rho_e_, i=1, 5)], [(i, i=0, 4)])
   call Q%print()

   call conservative%set_eos(eos)
   call conservative%push(Q)

   call domain1%create([5], dof=1, stencil_width=1)
   call u_vec%create(domain1)
   call v_vec%create(domain1)
   call w_vec%create(domain1)
   call rho_vec%create(domain1)
   call T_vec%create(domain1)

   call primitive%set_eos(eos)
   call primitive%push(rho_vec, u_vec, v_vec, w_vec, T_vec)

   !!! p to c
   call transform_field_c2p(conservative, primitive)
   ! call transform_field(primitive, conservative)
   call primitive%rho%print()
   call primitive%u%print()
   call primitive%v%print()
   call primitive%w%print()
   call primitive%T%print()

   call transform_field_p2c(primitive, conservative)
   call conservative%q%print()

   call conservative%destroy()
   call u_vec%destroy()
   call v_vec%destroy()
   call w_vec%destroy()
   call rho_vec%destroy()
   call T_vec%destroy()
   call Q%destroy()
   call domain%destroy(); call domain1%destroy()
   call eos%destroy()
   call petsc_finalize()

end program ex_field_conservative
